if(!require(ggplot2))install.packages("ggplot2"); library(ggplot2)
## Loading required package: ggplot2
if(!require(ggmap)) install.packages("ggmap"); library(ggmap)
## Loading required package: ggmap
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
# ggmap(get_map(location='south korea', zoom=7))
#error 발생 : google maps API terms of Service가 변경됨
if(!requireNamespace("devtools")) install.packages("devtools")
## Loading required namespace: devtools
devtools::install_github("dkahle/ggmap", ref = "tidyup")
## Skipping install of 'ggmap' from a github remote, the SHA1 (2d756e5e) has not changed since last install.
## Use `force = TRUE` to force installation
# 좌, 우, 위, 아래의 값을 지정해 주어야 지도 그림이 그려짐
kor <- c(left = 124, bottom = 33, right = 132, top = 39)
map <- get_stamenmap(kor, zoom = 6, maptype = "toner-lite")
## Source : http://tile.stamen.com/toner-lite/6/54/24.png
## Source : http://tile.stamen.com/toner-lite/6/55/24.png
## Source : http://tile.stamen.com/toner-lite/6/54/25.png
## Source : http://tile.stamen.com/toner-lite/6/55/25.png
ggmap(map)

map <- get_stamenmap(kor,
zoom = 6, maptype = "watercolor")
## Source : http://tile.stamen.com/watercolor/6/54/24.jpg
## Source : http://tile.stamen.com/watercolor/6/55/24.jpg
## Source : http://tile.stamen.com/watercolor/6/54/25.jpg
## Source : http://tile.stamen.com/watercolor/6/55/25.jpg
ggmap(map)

if(!require(dplyr)) install.packages("dplyr");library(dplyr)
## Loading required package: dplyr
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
if(!require(forcats)) install.packages("forcats");library(forcats)
## Loading required package: forcats
setwd("C:/Users/USER/Dropbox/R-work/vis/day2/extra")
load(file='crime.rda')
head(crime)
## time date hour premise offense beat
## 82729 2010-01-01 15:00:00 1/1/2010 0 18A murder 15E30
## 82730 2010-01-01 15:00:00 1/1/2010 0 13R robbery 13D10
## 82731 2010-01-01 15:00:00 1/1/2010 0 20R aggravated assault 16E20
## 82732 2010-01-01 15:00:00 1/1/2010 0 20R aggravated assault 2A30
## 82733 2010-01-01 15:00:00 1/1/2010 0 20A aggravated assault 14D20
## 82734 2010-01-01 15:00:00 1/1/2010 0 20R burglary 18F60
## block street type suffix number month day
## 82729 9600-9699 marlive ln - 1 january friday
## 82730 4700-4799 telephone rd - 1 january friday
## 82731 5000-5099 wickview ln - 1 january friday
## 82732 1000-1099 ashland st - 1 january friday
## 82733 8300-8399 canyon - 1 january friday
## 82734 9300-9399 rowan ln - 1 january friday
## location address lon lat
## 82729 apartment parking lot 9650 marlive ln -95.43739 29.67790
## 82730 road / street / sidewalk 4750 telephone rd -95.29888 29.69171
## 82731 residence / house 5050 wickview ln -95.45586 29.59922
## 82732 residence / house 1050 ashland st -95.40334 29.79024
## 82733 apartment 8350 canyon -95.37791 29.67063
## 82734 residence / house 9350 rowan ln -95.54830 29.70223
table(crime$offense)
##
## aggravated assault auto theft burglary
## 7177 7946 17802
## murder rape robbery
## 157 378 6298
## theft
## 46556
# define helper
`%notin%` <- function(lhs, rhs) !(lhs %in% rhs)
# reduce crime to violent crimes in downtown houston
violent_crimes <- crime %>%
filter(
offense %notin% c("auto theft", "theft", "burglary"),
-95.39681 <= lon & lon <= -95.34188,
29.73631 <= lat & lat <= 29.78400
) %>%
mutate(
offense = fct_drop(offense), # refactor
offense = fct_relevel(offense,
c("robbery", "aggravated assault", "rape", "murder") #factor order
)
)
head(violent_crimes)
## time date hour premise offense beat
## 1 2010-01-01 20:00:00 1/1/2010 5 20R aggravated assault 2A10
## 2 2010-01-01 21:00:00 1/1/2010 6 20R rape 10H10
## 3 2010-01-03 09:00:00 1/2/2010 18 18T robbery 1A20
## 4 2010-01-04 05:00:00 1/3/2010 14 210 robbery 6B10
## 5 2010-01-05 19:00:00 1/5/2010 4 20A aggravated assault 2A30
## 6 2010-01-06 05:00:00 1/5/2010 14 070 aggravated assault 10H40
## block street type suffix number month day
## 1 2000-2099 keene - 1 january friday
## 2 2300-2399 sperber ln - 1 january friday
## 3 400-499 gray W 1 january saturday
## 4 5500-5599 north fwy ser E 1 january sunday
## 5 2200-2299 white oak dr - 1 january tuesday
## 6 1000-1099 alabama st - 1 january tuesday
## location address lon lat
## 1 residence / house 2050 keene -95.36348 29.77802
## 2 residence / house 2350 sperber ln -95.34817 29.75505
## 3 strip business center parking lot 450 gray -95.37596 29.75143
## 4 restaurant / cafeteria 5550 north fwy ser -95.36327 29.76328
## 5 apartment 2250 white oak dr -95.38217 29.78002
## 6 convenience store 1050 alabama st -95.37930 29.73745
table(violent_crimes$offense)
##
## robbery aggravated assault rape
## 312 362 22
## murder
## 5
# use qmplot to make a scatterplot on a map
qmplot(lon, lat, data = violent_crimes, maptype = "toner-lite", color = I("red"))
## Using zoom = 14...
## Source : http://tile.stamen.com/terrain/14/3850/6770.png
## Source : http://tile.stamen.com/terrain/14/3851/6770.png
## Source : http://tile.stamen.com/terrain/14/3852/6770.png
## Source : http://tile.stamen.com/terrain/14/3853/6770.png
## Source : http://tile.stamen.com/terrain/14/3850/6771.png
## Source : http://tile.stamen.com/terrain/14/3851/6771.png
## Source : http://tile.stamen.com/terrain/14/3852/6771.png
## Source : http://tile.stamen.com/terrain/14/3853/6771.png
## Source : http://tile.stamen.com/terrain/14/3850/6772.png
## Source : http://tile.stamen.com/terrain/14/3851/6772.png
## Source : http://tile.stamen.com/terrain/14/3852/6772.png
## Source : http://tile.stamen.com/terrain/14/3853/6772.png
## Source : http://tile.stamen.com/terrain/14/3850/6773.png
## Source : http://tile.stamen.com/terrain/14/3851/6773.png
## Source : http://tile.stamen.com/terrain/14/3852/6773.png
## Source : http://tile.stamen.com/terrain/14/3853/6773.png

qmplot(lon, lat, data = violent_crimes, maptype = "toner-lite", geom = "density2d", color = I("red"))
## Using zoom = 14...

robberies <- violent_crimes %>% filter(offense == "robbery")
qmplot(lon, lat, data = violent_crimes, geom = "blank",
zoom = 15, maptype = "toner-background", darken = .7, legend = "topleft") +
stat_density_2d(aes(fill = ..level..), geom = "polygon", alpha = .3, color = NA)+
scale_fill_gradient2("Robbery\nPropensity", low = "white", mid = "yellow", high = "red", midpoint = 650)
## 49 tiles needed, this may take a while (try a smaller zoom).
## Source : http://tile.stamen.com/terrain/15/7700/13541.png
## Source : http://tile.stamen.com/terrain/15/7701/13541.png
## Source : http://tile.stamen.com/terrain/15/7702/13541.png
## Source : http://tile.stamen.com/terrain/15/7703/13541.png
## Source : http://tile.stamen.com/terrain/15/7704/13541.png
## Source : http://tile.stamen.com/terrain/15/7705/13541.png
## Source : http://tile.stamen.com/terrain/15/7706/13541.png
## Source : http://tile.stamen.com/terrain/15/7700/13542.png
## Source : http://tile.stamen.com/terrain/15/7701/13542.png
## Source : http://tile.stamen.com/terrain/15/7702/13542.png
## Source : http://tile.stamen.com/terrain/15/7703/13542.png
## Source : http://tile.stamen.com/terrain/15/7704/13542.png
## Source : http://tile.stamen.com/terrain/15/7705/13542.png
## Source : http://tile.stamen.com/terrain/15/7706/13542.png
## Source : http://tile.stamen.com/terrain/15/7700/13543.png
## Source : http://tile.stamen.com/terrain/15/7701/13543.png
## Source : http://tile.stamen.com/terrain/15/7702/13543.png
## Source : http://tile.stamen.com/terrain/15/7703/13543.png
## Source : http://tile.stamen.com/terrain/15/7704/13543.png
## Source : http://tile.stamen.com/terrain/15/7705/13543.png
## Source : http://tile.stamen.com/terrain/15/7706/13543.png
## Source : http://tile.stamen.com/terrain/15/7700/13544.png
## Source : http://tile.stamen.com/terrain/15/7701/13544.png
## Source : http://tile.stamen.com/terrain/15/7702/13544.png
## Source : http://tile.stamen.com/terrain/15/7703/13544.png
## Source : http://tile.stamen.com/terrain/15/7704/13544.png
## Source : http://tile.stamen.com/terrain/15/7705/13544.png
## Source : http://tile.stamen.com/terrain/15/7706/13544.png
## Source : http://tile.stamen.com/terrain/15/7700/13545.png
## Source : http://tile.stamen.com/terrain/15/7701/13545.png
## Source : http://tile.stamen.com/terrain/15/7702/13545.png
## Source : http://tile.stamen.com/terrain/15/7703/13545.png
## Source : http://tile.stamen.com/terrain/15/7704/13545.png
## Source : http://tile.stamen.com/terrain/15/7705/13545.png
## Source : http://tile.stamen.com/terrain/15/7706/13545.png
## Source : http://tile.stamen.com/terrain/15/7700/13546.png
## Source : http://tile.stamen.com/terrain/15/7701/13546.png
## Source : http://tile.stamen.com/terrain/15/7702/13546.png
## Source : http://tile.stamen.com/terrain/15/7703/13546.png
## Source : http://tile.stamen.com/terrain/15/7704/13546.png
## Source : http://tile.stamen.com/terrain/15/7705/13546.png
## Source : http://tile.stamen.com/terrain/15/7706/13546.png
## Source : http://tile.stamen.com/terrain/15/7700/13547.png
## Source : http://tile.stamen.com/terrain/15/7701/13547.png
## Source : http://tile.stamen.com/terrain/15/7702/13547.png
## Source : http://tile.stamen.com/terrain/15/7703/13547.png
## Source : http://tile.stamen.com/terrain/15/7704/13547.png
## Source : http://tile.stamen.com/terrain/15/7705/13547.png
## Source : http://tile.stamen.com/terrain/15/7706/13547.png

qmplot(lon, lat, data = violent_crimes, maptype = "toner-background", color = offense) +
facet_wrap(~ offense)
## Using zoom = 14...

#유럽지도 가져오기
europe <- c(left = -12, bottom = 35, right = 30, top = 63)
get_stamenmap(europe, zoom = 5) %>% ggmap()
## Source : http://tile.stamen.com/terrain/5/14/8.png
## Source : http://tile.stamen.com/terrain/5/15/8.png
## Source : http://tile.stamen.com/terrain/5/16/8.png
## Source : http://tile.stamen.com/terrain/5/17/8.png
## Source : http://tile.stamen.com/terrain/5/18/8.png
## Source : http://tile.stamen.com/terrain/5/14/9.png
## Source : http://tile.stamen.com/terrain/5/15/9.png
## Source : http://tile.stamen.com/terrain/5/16/9.png
## Source : http://tile.stamen.com/terrain/5/17/9.png
## Source : http://tile.stamen.com/terrain/5/18/9.png
## Source : http://tile.stamen.com/terrain/5/14/10.png
## Source : http://tile.stamen.com/terrain/5/15/10.png
## Source : http://tile.stamen.com/terrain/5/16/10.png
## Source : http://tile.stamen.com/terrain/5/17/10.png
## Source : http://tile.stamen.com/terrain/5/18/10.png
## Source : http://tile.stamen.com/terrain/5/14/11.png
## Source : http://tile.stamen.com/terrain/5/15/11.png
## Source : http://tile.stamen.com/terrain/5/16/11.png
## Source : http://tile.stamen.com/terrain/5/17/11.png
## Source : http://tile.stamen.com/terrain/5/18/11.png
## Source : http://tile.stamen.com/terrain/5/14/12.png
## Source : http://tile.stamen.com/terrain/5/15/12.png
## Source : http://tile.stamen.com/terrain/5/16/12.png
## Source : http://tile.stamen.com/terrain/5/17/12.png
## Source : http://tile.stamen.com/terrain/5/18/12.png

kor <- c(left = 124, bottom = 33, right = 132, top = 39)
map <- get_stamenmap(kor, zoom = 7, maptype = "toner-lite")
## Source : http://tile.stamen.com/toner-lite/7/108/48.png
## Source : http://tile.stamen.com/toner-lite/7/109/48.png
## Source : http://tile.stamen.com/toner-lite/7/110/48.png
## Source : http://tile.stamen.com/toner-lite/7/108/49.png
## Source : http://tile.stamen.com/toner-lite/7/109/49.png
## Source : http://tile.stamen.com/toner-lite/7/110/49.png
## Source : http://tile.stamen.com/toner-lite/7/108/50.png
## Source : http://tile.stamen.com/toner-lite/7/109/50.png
## Source : http://tile.stamen.com/toner-lite/7/110/50.png
## Source : http://tile.stamen.com/toner-lite/7/108/51.png
## Source : http://tile.stamen.com/toner-lite/7/109/51.png
## Source : http://tile.stamen.com/toner-lite/7/110/51.png
plot(map)

wifi <- read.csv('wifi.csv', header=T) #wifi.csv
head(wifi)
## company lat lon
## 1 KT 37.74417 128.9056
## 2 KT 37.72806 128.9543
## 3 KT 37.75710 128.8900
## 4 KT 37.74769 128.8840
## 5 KT 37.74866 128.9073
## 6 KT 37.74281 128.8827
ggmap(map) + geom_point(data=wifi, aes(x=lon, y=lat, color=company))

ggmap(map) + stat_density_2d(data=wifi, aes(x=lon, y=lat))

ggmap(map) +
stat_density_2d(data=wifi, aes(x=lon, y=lat, fill=..level.., alpha=..level..), geom='polygon', size=2, bins=30)

#일단 aes 안에 들어 있는 fill은 문자 그대로 색깔로 채우라는 뜻입니다.
# ..level..은 레벨(level)이 높을수록, 예를 들어 기압이 높거나 고도가 높을수록
# 더 진한 색깔을 칠하라는 뜻입니다.
# alpha는 투명도를 나타냅니다.
# 역시 같은 원리로 레벨이 높으면 불투명하게(색이 더 잘 드러나게)
# 칠하고 낮을 때는 투명하게(희미하게) 칠하라는 의미입니다
# geom='polygon'에서 polygon은 다각형이라는 뜻입니다.
# 기본은 위에서 본 것처럼 선(線)으로 돼 있는데 그것 말고 도형으로 그리라고 명령을 준 겁니다.
# 선을 기준으로 size는 선 굵기, bins는 선 간격이라는 뜻입니다.
# 이름이 이렇게 붙은 건 원래 점을 이은 선을 기준으로 삼고 있는 까닭입니다.
p <- ggmap(map) + stat_density_2d(data=wifi, aes(x=lon, y=lat, fill=..level.., alpha=..level..), geom='polygon', size=7, bins=28)
p

p<-p + scale_fill_gradient(low='yellow', high='red', guide=F)
p

p<-p + scale_alpha(range=c(0.02, 0.8), guide=F)
p

airport <- read.csv('airport.csv', header=T) #airport.csv
route <- read.csv('route.csv', header=T) #route.csv
head(airport)
## airport iata lon lat
## 1 강릉 KAG 128.944 37.7536
## 2 광주 KWJ 126.809 35.1264
## 3 군산 KUV 126.616 35.9038
## 4 김포 GMP 126.791 37.5583
## 5 대구 TAE 128.659 35.8941
## 6 목포 MPK 126.380 34.7589
head(route)
## id airport lon lat
## 1 1 CJJ 127.499 36.7166
## 2 7 CJJ 127.499 36.7166
## 3 45 CJJ 127.499 36.7166
## 4 77 CJJ 127.499 36.7166
## 5 2 CJJ 127.499 36.7166
## 6 8 CJJ 127.499 36.7166
ggmap(map) + geom_point(data=airport, aes(x=lon, y=lat, size=3))

p <- ggmap(map) + geom_point(data=airport, aes(x=lon, y=lat, size=3))
p<-p + geom_line(data=route, aes(x=lon, y=lat, group=id))
p

CJU_route<-route[route$airport=="CJU","id"]
route1 <- route[route$id %in% CJU_route,]
p + geom_line(data=route1, aes(x=lon, y=lat, group=id, size=2), col='gold')

# Q. 기상관측소 데이터를 찾아서 지도에 그려보세요.
# stn<- read.csv(choose.files())
if(!require(raster))install.packages("raster"); library(raster)
## Loading required package: raster
## Loading required package: sp
##
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
##
## select
korea <- getData('GADM', country='kor', level=2)
plot(korea)

ggplot(data=korea, aes(x=long, y=lat, group=group)) +
geom_polygon(fill='white', color='black')
## Regions defined for each Polygons

# 잘 보시면 예전에 못 보던 게 하나 보입니다.
# group : 선을 모아 면을 만들 때도 짝을 맞춰주는 속성
# 이 지도는 옛날 행정구역임
if(!require(rgdal))install.packages("rgdal"); library(rgdal)
## Loading required package: rgdal
## rgdal: version: 1.4-4, (SVN revision 833)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/USER/Documents/R/win-library/3.6/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/USER/Documents/R/win-library/3.6/rgdal/proj
## Linking to sp version: 1.3-1
if(!require(maptools))install.packages("maptools"); library(maptools)
## Loading required package: maptools
## Checking rgeos availability: FALSE
## Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,
## which has a restricted licence. It is disabled by default;
## to enable gpclib, type gpclibPermit()
korea <- readOGR('SIG_201703/TL_SCCO_SIG.shp')#'TL_SCCO_SIG.shp'
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\USER\Dropbox\R-work\vis\day2\extra\SIG_201703\TL_SCCO_SIG.shp", layer: "TL_SCCO_SIG"
## with 250 features
## It has 3 fields
ggplot() + geom_polygon(data=korea, aes(x=long, y=lat, group=group), fill='white', color='black')
## Regions defined for each Polygons

#19대 대통령 선거 시군구별 득표율
result <- read.csv('result.csv', header=T) #'result.csv'
head(result)
## Sido Sigun id moon hong ahn yoo shim etc X1st X1st_C
## 1 서울 종로구 11110 0.416 0.218 0.218 0.073 0.070 0.005 0.416 Moon
## 2 서울 중구 11140 0.412 0.217 0.235 0.071 0.060 0.005 0.412 Moon
## 3 서울 용산구 11170 0.393 0.239 0.217 0.080 0.066 0.004 0.393 Moon
## 4 서울 성동구 11200 0.428 0.200 0.226 0.078 0.064 0.004 0.428 Moon
## 5 서울 광진구 11215 0.441 0.194 0.221 0.072 0.069 0.004 0.441 Moon
## 6 서울 동대문구 11230 0.421 0.219 0.227 0.064 0.064 0.005 0.421 Moon
## X2nd X2nd_C gap final
## 1 0.218 Hong 0.198 0.416
## 2 0.235 Ahn 0.178 0.412
## 3 0.239 Hong 0.155 0.393
## 4 0.226 Ahn 0.203 0.428
## 5 0.221 Ahn 0.220 0.441
## 6 0.227 Ahn 0.194 0.421
head(korea)
## SIG_CD SIG_ENG_NM SIG_KOR_NM
## 0 11110 Jongno-gu 종로구
## 1 11140 Jung-gu 중구
## 2 11170 Yongsan-gu 용산구
## 3 11200 Seongdong-gu 성동구
## 4 11215 Gwangjin-gu 광진구
## 5 11230 Dongdaemun-gu 동대문구
korea1<-merge(korea, result, by.x="SIG_CD", by.y="id")
head(korea1)
## SIG_CD SIG_ENG_NM SIG_KOR_NM Sido Sigun moon hong ahn yoo
## 1 11110 Jongno-gu 종로구 서울 종로구 0.416 0.218 0.218 0.073
## 2 11140 Jung-gu 중구 서울 중구 0.412 0.217 0.235 0.071
## 3 11170 Yongsan-gu 용산구 서울 용산구 0.393 0.239 0.217 0.080
## 4 11200 Seongdong-gu 성동구 서울 성동구 0.428 0.200 0.226 0.078
## 5 11215 Gwangjin-gu 광진구 서울 광진구 0.441 0.194 0.221 0.072
## 6 11230 Dongdaemun-gu 동대문구 서울 동대문구 0.421 0.219 0.227 0.064
## shim etc X1st X1st_C X2nd X2nd_C gap final
## 1 0.070 0.005 0.416 Moon 0.218 Hong 0.198 0.416
## 2 0.060 0.005 0.412 Moon 0.235 Ahn 0.178 0.412
## 3 0.066 0.004 0.393 Moon 0.239 Hong 0.155 0.393
## 4 0.064 0.004 0.428 Moon 0.226 Ahn 0.203 0.428
## 5 0.069 0.004 0.441 Moon 0.221 Ahn 0.220 0.441
## 6 0.064 0.005 0.421 Moon 0.227 Ahn 0.194 0.421
SIG_CD<-paste(unique(korea1$SIG_CD))
f.korea <- fortify(korea) #maptools package
## Regions defined for each Polygons
f.korea$SIG_CD <-NA
for (i in 1:250){
f.korea[f.korea$id==(i-1),"SIG_CD"]=SIG_CD[i]
}
korea2 <- merge(f.korea, result, by.x="SIG_CD", by.y="id")
p<-ggplot() + geom_polygon(data=korea2, aes(x=long, y=lat, group=group, fill=moon))
p + scale_fill_gradient(low='white', high='#004ea2')

if(!require(viridis))install.packages("viridis");library(viridis)
## Loading required package: viridis
## Loading required package: viridisLite
p<-ggplot() + geom_polygon(data=korea2, aes(x=long, y=lat, group=group, fill=moon))
p + scale_fill_viridis()

# 기존에 있는 파레트를 이용해서 그려보자.
p + scale_fill_viridis(direction=-1) + theme_void() + guides(fill=F)

# 색을 반대로, 회색을 지우고, 범례도 지워보자
# 다른 후보들의 그림을 그려보세요.
# 문재인 후보와 홍준표 후보의 차이를 시각적으로 살펴보세요.